doseresNMA:

#!!! CHECK: knots <- quantile(min,max) OR quantile(0,max)???
net.structure <- function(data=data,
                          metareg=F,class.effect=F,
                          ref.lab='placebo',refclass.lab='placebo',
                          knot_probs=c(0.25,0.50,0.75)
){
  # 1. data: data.frame() that have the following
  # studyid
  # r: number of events
  # n: sample size
  # dose: the dose level
  # drug: drug names
  # cov (if metareg is TRUE)
  # class (if class.effect is TRUE)
  # 2. metareg: T/F to preform dose-effect NMR
  # 3. class: T/F to run dose-effect class model
  # 4. ref.lab: the referenced drug, default is placebo
  # 5. refclass.lab: the referenced class, default is placebo
  
  # load
  require('dplyr')
  require('rms')
  #** initial arguments
  ns <- length(unique(data$studyid))    # number of studies
  ndrugs <- length(unique(data$drug))   # number of drugs
  na <- as.numeric(table(data$studyid)) # number of arms per study
  max.na <- max(na)                     # maximum number of arms
  
  data$studyID <- as.numeric(as.factor(data$studyid)) # add numeric odered study id
  study_id <- unique(data$studyID)      # a vector of unique studyID
  
  data$drug_index <- as.numeric(factor(data$drug)) # drug numeric index
  data$drug <- tolower(data$drug) # translate drug names to lower case
  
  # the numeric index of the reference
  if(is.null(ref.lab)&is.null(refclass.lab)){
    ref.lab<- refclass.lab <-'placebo'
  }else{
    ref.lab<-ref.lab
    refclass.lab <- refclass.lab
  }
  ref.index <- unique(data$drug_index[data$drug==ref.lab])
  refclass.index <- unique(data$classF[data$class==refclass.lab])
  
  #** doses to RCS transformation
  
  # dataset without placebo to obtain RCS transformation for each non-placebo drug
  data_no_placebo <- data[data$drug!=ref.lab,]
  data_no_placebo$drug <- factor(data_no_placebo$drug)
  
  if(class.effect){ # for hramonized doses
    # FIND RCS
    max.dose <- max(data$dose,na.rm = T)
    min.dose <- min(data$dose,na.rm = T)
    knots <- quantile(min.dose:max.dose,probs = knot_probs)
    rcsdose <-  rcspline.eval(data$dose,knots = knots,inclx = TRUE)
    
    # ADD RCS
    data$dose2 <- rcsdose[,2]
    
    # FIND RCS for absolute predictions
    new.dose <-  1:max.dose
    f.new.dose <-  rcspline.eval(new.dose,knots = knots,inclx = TRUE)[,2]
    nd.new <- round(max.dose)
    
  }else{ # for unhramonized doses
    # FIND rcs
    rcsdose_drug <- data_no_placebo %>%
      group_by(drug) %>%
      group_map(~dose_to_rcs(.x$dose)) # the list order as the order of levels(data_no_placebo)
    rcsdose <- do.call(rbind,rcsdose_drug)
    rownames(rcsdose) <- rep(levels(data_no_placebo$drug),
                             times=as.numeric(table(data_no_placebo$drug)))
    # ADD rcsdose: dose2
    data$dose2 <- data$dose
    for (i in 1:(ndrugs-1)){
      data$dose2[data$drug==levels(data_no_placebo$drug)[i]] <- rcsdose[rownames(rcsdose)==levels(data_no_placebo$drug)[i],2]
    }
    # FIND RCS for absolute predictions
    max.dose <- round(max(data$dose,na.rm = T))
    new.dose <- f.new.dose <- matrix(NA,ndrugs,max.dose)
    nd.new <-  rep(0,ndrugs)
    
    for (i in unique(data_no_placebo$drug_index)){ # loop through non-placebo drugs
      max.dose <- round(max(data[data$drug_index==i,]$dose,na.rm = T))
      min.dose <- round(min(data[data$drug_index==i,]$dose,na.rm = T))
      new.dose[i,1:max.dose] <-  1:max.dose 
      knots <- quantile(min.dose:max.dose,probs = knot_probs)
      f.new.dose[i,1:max.dose] <-  rcspline.eval(new.dose[i,1:max.dose],knots = knots,inclx = TRUE)[,2]
      nd.new[i] <- length(new.dose[i,1:max.dose])
    }
    
  }
  # ORDER the doses on each study: start from zero
  data_withRCS <-data%>%
    group_by(studyID)%>%
    arrange(dose,.by_group=TRUE)
  
  #** Column to  matrix with ns in rows and na in columns(max.na is specified, the additionals are given NA)
  rmat <- col_to_mat(data_withRCS,data_withRCS$r)
  
  nmat <- col_to_mat(data_withRCS,data_withRCS$n)
  
  dosemat <- col_to_mat(data_withRCS,data_withRCS$dose)
  
  dose2mat <- col_to_mat(data_withRCS,data_withRCS$dose2)
  
  tmat <- col_to_mat(data_withRCS,data_withRCS$drug_index)
  
  #** For inconsistency model: indices of direct comparisons
  comp <- direct.comp.index(data_withRCS)
  t1 <- comp[,'t1']
  t2 <- comp[,'t2']
  ncomp <- nrow(comp)
  
  #** RETURN: list to jags models
  jagsdata <- list(ns=ns, na=na,ndrugs=ndrugs,
                   r=rmat, n=nmat,dose1=dosemat,t=tmat,dose2=dose2mat,
                   ref.index=ref.index,
                   rr=rmat,nn=nmat,new.dose=new.dose,f.new.dose=f.new.dose,nd.new=nd.new,
                   t1=t1,t2=t2,ncomp=ncomp
  )
  #** Add to jagsdata for metareg or class effect parts 
  if(metareg){
    if(class.effect){
      classmat <-col_to_mat(data_withRCS,data_withRCS$classF)
      nc <- length(unique(classF))
      covmat <-col_to_mat(data_withRCS,data_withRCS$cov)
      add <- list(class.effect=classmat, refclass.index=refclass.index, nc=nc,
                  cov=colMeans(t(covmat),na.rm=T)
      )
      jagsdata <- c(jagsdata,add)
    }else{
      covmat <-col_to_mat(data_withRCS,data_withRCS$cov)
      add <- list(cov=colMeans(t(covmat),na.rm=T)
      )
      jagsdata <- c(jagsdata,add)
    }
  }else{
    if(class.effect){
      classmat <-col_to_mat(data_withRCS,data_withRCS$classF)
      nc <- length(unique(data_withRCS$classF))
      add <- list(class=classmat, refclass.index=refclass.index, nc=nc)
      jagsdata <- c(jagsdata,add)
    }else{
      jagsdata
    }
  }
  
  return(jagsdata)
}


#!!! think about it and ask Georgia: the quantile should be computed from min OR 0 to max dose
# functions to be used in net.structure

# 1. take the drug-dose and find its RCS transformation at
dose_to_rcs <- function(dose.per.drug,knot_probs=c(0.25,0.50,0.75)){
  require('rms')
  max.dose <- max(dose.per.drug)
  min.dose <- min(dose.per.drug)
  knots <- quantile(0:max.dose,probs = knot_probs)
  rcs.dose.per.drug <-  rcspline.eval(dose.per.drug,knots = knots,inclx = TRUE)
  return(rcs.dose.per.drug)
}

# 2. convert the data column to matrix ( row is study and column is arm (dose-level) )

col_to_mat <- function(data,var){
  ns <-length(unique(data$studyid))
  na <- as.numeric(table(data$studyid)) # number of arms per study
  max.na <- max(na)
  data$studyID <- as.numeric(as.factor(data$studyid))         # transform studyid to ordered numeric values
  study_id <- unique(data$studyID)
  varmat <- matrix(NA,ns,max.na)
  for (i in 1:ns) {
    varmat[i,1:as.numeric(table(data$studyID)[i])] <- var[data$studyID == study_id[i]]
  }
  return(varmat)
}

# 3. find the indices of the direct head-to-head comparisons
direct.comp.index <- function(data)
{
  require(dplyr)
  data <- dplyr::arrange(data, data$studyid, data$dose)
  t1 <- vector()
  t2 <- vector()
  for (i in seq_along(unique(data[["studyid"]]))) {
    subset <- subset(data, studyid==unique(data[["studyid"]])[i])
    for (k in 2:nrow(subset)) {
      t1 <- append(t1, subset[["drug"]][1])
      t2 <- append(t2, subset[["drug"]][k])
      if (is.na(subset[["drug"]][k])) {
        stop()
      }
    }
  }
  
  comparisons <- data.frame(t1 = t1, t2 = t2)
  comparisons <- comparisons %>% dplyr::group_by(t1, t2) %>%
    dplyr::mutate(nr = dplyr::n())
  comparisons <- unique(comparisons)
  comparisons <- dplyr::arrange(comparisons, t1, t2)
  row_name = comparisons$row_name
  comparisons %<>% select(-row_name) %>% as.matrix
  rownames(comparisons) = row_name
  return(comparisons)
}